﻿using System;

namespace OpenLS.Spreadsheet.StandardAddIn
{
    
   
    internal static class Beta {
        private static readonly double[] Y = {
                                       0.0021695375159141994,
                                       0.011413521097787704, 0.027972308950302116, 0.051727015600492421,
                                       0.082502225484340941, 0.12007019910960293, 0.16415283300752470,
                                       0.21442376986779355, 0.27051082840644336, 0.33199876341447887,
                                       0.39843234186401943, 0.46931971407375483, 0.54413605556657973,
                                       0.62232745288031077, 0.70331500465597174, 0.78649910768313447,
                                       0.87126389619061517, 0.95698180152629142
                                   };

        private static readonly double[] W = {
                                       0.0055657196642445571,
                                       0.012915947284065419, 0.020181515297735382, 0.027298621498568734,
                                       0.034213810770299537, 0.040875750923643261, 0.047235083490265582,
                                       0.053244713977759692, 0.058860144245324798, 0.064039797355015485,
                                       0.068745323835736408, 0.072941885005653087, 0.076598410645870640,
                                       0.079687828912071670, 0.082187266704339706, 0.084078218979661945,
                                       0.085346685739338721, 0.085983275670394821
                                   };

        private static readonly double[] Cof = {
                                                   57.1562356658629235, -59.5979603554754912,
                                                   14.1360979747417471, -0.491913816097620199, .339946499848118887e-4,
                                                   .465236289270485756e-4, -.983744753048795646e-4, .158088703224912494e-3
                                                   ,
                                                   -.210264441724104883e-3, .217439618115212643e-3,
                                                   -.164318106536763890e-3,
                                                   .844182239838527433e-4, -.261908384015814087e-4, .368991826595316234e-5
                                               };

        private const int Switch = 3000;

        private const double Fpmin = double.MinValue/double.Epsilon;


        public static double BetaI(double a, double b, double x)
        {
            if (a <= 0.0 ) throw new ArgumentException("a");
            if (b <= 0.0) throw new ArgumentException("b");
            if (x < 0.0 || x > 1.0) throw new ArgumentException("x");
            if (x == 0.0 || x == 1.0) return x;
            if (a > Switch && b > Switch) return betaiapprox(a, b, x);
            double bt = Math.Exp(gammln(a + b) - gammln(a) - gammln(b) + a*Math.Log(x) + b*Math.Log(1.0 - x));
            if (x < (a + 1.0)/(a + b + 2.0)) return bt*betacf(a, b, x)/a;
            return 1.0 - bt*betacf(b, a, 1.0 - x)/b;
        }

        private static double betacf(double a, double b, double x)
        {
            int m;
            double qab = a + b;
            double qap = a + 1.0;
            double qam = a - 1.0;
            double c = 1.0;
            double d = 1.0 - qab*x/qap;
            if (Math.Abs(d) < Fpmin) d = Fpmin;
            d = 1.0/d;
            double h = d;
            for (m = 1; m < 10000; m++)
            {
                int m2 = 2*m;
                double aa = m*(b - m)*x/((qam + m2)*(a + m2));
                d = 1.0 + aa*d;
                if (Math.Abs(d) < Fpmin) d = Fpmin;
                c = 1.0 + aa/c;
                if (Math.Abs(c) < Fpmin) c = Fpmin;
                d = 1.0/d;
                h *= d*c;
                aa = -(a + m)*(qab + m)*x/((a + m2)*(qap + m2));
                d = 1.0 + aa*d;
                if (Math.Abs(d) < Fpmin) d = Fpmin;
                c = 1.0 + aa/c;
                if (Math.Abs(c) < Fpmin) c = Fpmin;
                d = 1.0/d;
                double del = d*c;
                h *= del;
                if (Math.Abs(del - 1.0) <= double.Epsilon) break;
            }
            return h;
        }


        private static double betaiapprox(double a, double b, double x)
        {
            int j;
            double xu;
            double a1 = a - 1.0, b1 = b - 1.0, mu = a / (a + b);
            double lnmu = Math.Log(mu), lnmuc = Math.Log(1.0 - mu);
            double d = a + b;
            double t = Math.Sqrt(a*b/(d*d*(a + b + 1.0)));
            if (x > a/(a + b))
            {
                if (x >= 1.0) return 1.0;
                xu = Math.Min(1.0, Math.Max(mu + 10.0*t, x + 5.0*t));
            }
            else
            {
                if (x <= 0.0) return 0.0;
                xu = Math.Max(0.0, Math.Min(mu - 10.0*t, x - 5.0*t));
            }
            double sum = 0;
            for (j = 0; j < 18; j++)
            {
                t = x + (xu - x)*Y[j];
                sum += W[j]*Math.Exp(a1*(Math.Log(t) - lnmu) + b1*(Math.Log(1 - t) - lnmuc));
            }
            double ans = sum*(xu - x)*Math.Exp(a1*lnmu - gammln(a) + b1*lnmuc - gammln(b) + gammln(a + b));
            return ans > 0.0 ? 1.0 - ans : -ans;
        }

        private static double gammln(double xx)
        {
            int j;
            double x;
            if (xx <= 0) throw new ArgumentException("xx");
            double y = x = xx;
            double tmp = x + 5.24218750000000000;
            tmp = (x + 0.5)*Math.Log(tmp) - tmp;
            double ser = 0.999999999999997092;
            for (j = 0; j < 14; j++) ser += Cof[j]/++y;
            return tmp + Math.Log(2.5066282746310005*ser/x);
        }
    }
}